capture program drop btstrp
program define btstrp, rclass
	
	tempvar propFxbt propfxbt polldayhatbt polldaybt Fxpeval
	egen `polldaybt'=sum(dummypollday)

	gb2lfit_mod pmconc, censvar(cens) from(intv)
	local esta = e(ba)
	local estb = e(bb)
	local estp = e(bp)
	local estq = e(bq)
	local atx=cff1
	
	gen `Fxpeval'=ibeta(`estp',`estq', (`atx'/`estb')^`esta'/(1+(`atx'/`estb')^`esta'))
	sum pmconc if pmconc<=float(cff1)
	gen `propFxbt'=r(N)/_N-`Fxpeval'  // total proportion of manipulation up to the first cutoff
	sum `propFxbt'
	return scalar propFxbtstrp=r(mean)
	return scalar estabtstrp=`esta'
	return scalar estbbtstrp=`estb'
	return scalar estpbtstrp=`estp'
	return scalar estqbtstrp=`estq'
	end

scalar lftcff1=0.135   
scalar rtcff1=0.18
scalar cff1=0.15

forv yearval=2001/2010{
use "./data/pm10data.dta", clear
di idval
di `yearval'
qui keep if id==idval
qui keep if year==`yearval'
di _N
//qui if _N>0{  // drop city/year if there are too many missing values

gen cens=0
replace cens=1 if pmconc>float(lftcff1) & pmconc <=float(cff1)  // if cens=1, these data are possibly manipulated so we censor them in estimation.
sum pmconc if pmconc<=float(cff1)
gen Fxnp=r(N)/_N  // Fx-np
sum pollday
if r(mean)>0{
DCdensity pmconc, breakpoint(0.15001) generate(Xj Yj r0 fhat se_fhat) nograph
drop Xj Yj r0 fhat se_fhat
gen mc=r(theta)
scalar h1=r(bandwidth)
}
else{
qui gen mc=.
}

gb2lfit_mod pmconc, cdf(pmparacdf) censvar(cens)

// the gb2lfit_mod codes only allow for manipulation of one direction

scalar lna0=ln(e(ba))
scalar lnb0=ln(e(bb))
scalar lnp0=ln(e(bp))
scalar lnq0=ln(e(bq))
matrix intv = (lna0, lnb0, lnp0, lnq0)
matrix colnames intv = ln_a:_cons ln_b:_cons ln_p:_cons ln_q:_cons

qui gen esta=e(ba)
qui gen estb=e(bb)
qui gen estp=e(bp)
qui gen estq=e(bq)
local ahat=e(ba)
local bhat=e(bb)
local phat=e(bp)
local qhat=e(bq)
local atx=cff1
gen Fxp=ibeta(`phat',`qhat', (`atx'/`bhat')^`ahat'/(1+(`atx'/`bhat')^`ahat'))

qui gen propFx=(Fxnp-Fxp)/Fxnp
qui gen polldayhat=(1-Fxp)*totday

bootstrap propFxbtstrp=r(propFxbtstrp)  ///
estabtstrp=r(estabtstrp) estbbtstrp=r(estbbtstrp) estpbtstrp=r(estpbtstrp) estqbtstrp=r(estqbtstrp), cluster(woy) ///
reps(200) saving("./data/bsreg.dta",replace): btstrp 
append using "./data/bsreg.dta"

egen estase=sd(estabtstrp)
egen estbse=sd(estbbtstrp)
egen estpse=sd(estpbtstrp)
egen estqse=sd(estqbtstrp)
egen propFxse=sd(propFxbtstrp) 
drop *btstrp cens 

keep in 1
local cityname=city[1] 

if `yearval'==2001{
save "./data/`cityname'_CMLEresults.dta", replace
} 
else {
qui merge 1:1 id year using "./data/`cityname'_CMLEresults.dta"  // merge dataset
drop _merge
qui save "./data/`cityname'_CMLEresults.dta", replace
}
}

scalar cv=invnorm(1-0.05/_N)
replace propFx=propFx*100  // make the unit be % change
replace propFxse=propFxse*100

gen propFxcileft=propFx-cv*propFxse
gen propFxciright=propFx+cv*propFxse

gen pctblueday=(totday-pollday)/totday
gen pctbluedayhat=(totday-polldayhat)/totday

// start plotting
local cityname=city[1]  
sort year
twoway (connected pctblueday year) (connected pctbluedayhat year), legend(order(1 "Reported" 2 "Estimated")) ///
xlabel(2002(2)2010) graphregion(color(white)) name(distest, replace) nodraw title("`cityname': % of Blue-sky Days")

twoway (connected propFx year) (line propFxcileft year, lcolor(gs1) lpattern(dash)) ///
(line propFxciright year,lcolor(gs1) lpattern(dash)), ///
legend(order(1 "Estimated %" 2 "90% Familywise CI")) xlabel(2002(2)2010) ylabel(`foo') graphregion(color(white)) ///
name(propFx, replace) nodraw title("% Manipulation Among Blue-sky Days")

